Scaling between Structural Relaxation and Particle Caging in a Model 
Colloidal Gel 
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In polymers melts and supercooled liquids, the glassy dynamics is characterized by the rattling of monomers or particles in 
the cage formed by their neighbors. Recently, a direct correlation in such systems, described by a universal scaling form, has 
been established between the rattling amplitude and the structural relaxation time. In this paper we analyze the glassy dynamics 
emerging from the formation of a persistent network in a model colloidal gel at very low density. The structural relaxation time 
of the gel network is compared with the mean squared displacement at short times, corresponding to the localization length 
associated to the presence of energetic bonds. Interestingly, we find that the same type of scaling as for the dense glassy systems 
holds. Our findings well elucidate the strong coupling between the cooperative rearrangements of the gel network and the single 
particle localization in the structure. Our results further indicate that the scaling captures indeed fundamental physical elements 
of glassy dynamics. 
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Understanding the extraordinary slow-down that accompanies 
systems with structural arrest like glass-forming systems^ 
and soft matterS^is a major scientific challenge. Crowding 
and caging effects play major roles in the glass transition (GT) 
of dense systems and lead to the strong localization of the par- 
ticles in the cage formed by their neighbors: this is typically 
apparent in the time dependence of the particle mean square 
displacement (MSD) as a plateau-like regime. 

The value of the MSD plateau (u 2 ) yields the amplitude of 
the rattling motion inside the cage and hence the correspond- 
ing localization length. 

Because of the extreme time-scale separation between the 
rattling motion (~ 10~ 12 s) and the structural relaxation (x a ~ 
10 2 s at GT), one would expect their complete independence. 

Nonetheless, there are several hints of the presence of cor- 
relations and several authors have investigated themE^H. i n 
particular, it has been recently shown that the structural re- 
laxation time T a and the rattling amplitude (m 2 ) 1 / 2 of several 
numerical models, including linear polymers, mixtures, pro- 
totypical glassformers like S1O2 and o-terphenyl (OTP), and 
one icosahedral glas sformer 25 " 27 !, can be related in a unique 
scaling form. 

Remarkably, the same resulting master curve well fits to the 
experimental data from van der Waals and associating liquids, 
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polymers, metallic glasses, ionic liquids and network glass- 
formers over many decades in timeESEI]. These results well 
elucidate how, in an impressively large class of systems, the 
glassy structural arrest corresponds to the onset of a strong 
coupling between the overall relaxation, characterized by co- 
operative and heterogeneous processes, and the average local- 
ization at the level of a single particle, strongly suggesting that 
this is a universal, fundamental feature of this type of dynam- 



Structural arrest and glassy dynamics can be observed also 
in very dilute particle suspensions when gelation occurs and 
displays significant hints of caging effects even at rather low 
volume fractionsEiKSl i n these cases, differently from the 
dense glass-forming systems just mentioned, most of the times 
the caging has been associated to the particle bonding, rather 
than to the role of excluded volume interactionsESHE] j n mese 
systems, particles get bonded into an interconnected network 
structure which is responsible for the onset of cooperative, 
slow dynamics and eventually structural arrest, i.e. gelation. 
It has been recently shown, in the numerical study of a model 
colloidal gel, that the gel network induces the same type of 
strong coupling in particle motion typically observed in dense 
glassy systems and that the glassy dynamics directly arises 
from the cooperative processes induced by the network^. A 
distinctive feature of this type of systems is the presence of 
different localization processes, over different length scales, 
leading to a somewhat more complex scenario for structural 
arrest. Nevertheless, particle caging is also observed, albeit 
much weaker than in dense systems, and again over time 
scales which are well separated from the ones of structural re- 
laxation. Intrigued by the similarities and differences in the 
glassy dynamics of dense glasses and low volume fraction 
gels, we have investigated the presence and nature of corre- 
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lations between structural relaxation and particle localization 
in a model colloidal gel. In spite of the deep differences in 
the caging mechanism and onset of slow dynamics, we find 
that the same universal scaling, already found in dense glassy 
systems, between the structural relaxation time x a and (u 2 ) 
holdsESEI]. j n Qur v j ew jjjjg fj nc [ m g points to the presence of 

a complex feedback between the overall slow structural relax- 
ation and the single particle localization in dilute gels. At the 
same time, it further supports the idea that this scaling form 
captures a fundamental, universal feature of glassy structural 
arrest. 




x (Configuration Space) 



Fig. 1 Two level systems where structural relaxation is achieved 
through a jump from one minimum to the other one overcoming an 
energy barrier AE (see text for details). 

The paper is organized as follows. In Sec. [2] we recall the 
basics to derive the universal scaling form subsequently dis- 
cussed, in Sec. [3] we briefly summarize the main features of 
the gel model as studied in Ref J 36 " 38 ! by molecular dynamics 
and of numerical simulations. In Sec. |4|we provide a sim- 
ple protocol to test the universal scaling^ 5 for the gel system 
studied here. Finally in Sec. [5] we draw the conclusions of our 
analysis. 



2 Universal Scaling Form 

On approaching the glass transition particles are longer and 
longer trapped into the cage formed by their neighbors. 
Caging phenomenon shows up as a plateau-like regime at 
short times in the MSD. The amplitude of the rattling motion 
(m 2 ) 1 / 2 during this caging regime, that occurs on very short 
time scales (e.g. picoseconds in molecular liquids), is the so- 
called Debye- Waller (DW) factor 16 l 2Q lwhich is directly related 
to the short-time elastic properties of the system^. 

The DW factor is an experimentally accessible quantity^! 
that can be also measured by using the incoherent interme- 




Fig. 2 Static structure factor at lowest temperatures investigated for 
all volume fractions. q m j„ and q max are are also pointed out where 
q m i„ is the minimum wave vector allowed by finite size of 
simulation box and q max is the wave vector corresponding to the 
maximum of S(q) at length scales comparable to particles diameter. 



diate scattering function (ISF), evaluating at the short times 
the height of the plateau which signals the cage effects (see 



Ref.^S). We note that as shown in^DW factor extracted from 
MSD (m 2 ) and one defined from ISF are equivalent. 

In spite of the fact that the DW factor is related to fast mo- 
tion of particles occurring on time scales much shorter than the 
ones typical of structural relaxation, many studies evidence in 
glass forming liquids a possible re lation between slow and fast 
degrees of feedonflIE™raEfl 

In order to express the correlation between DW factor and 
structural relaxation time in a functional form, a classical argu- 
ment estimating the height of the barrier between two poten- 
tial energy minima from the curvature around the minima can 
be used. For glassy systems Hall and Wolynes 13 applied this 
argument in their density functional theory where atomic mo- 
tion is restricted to cells, picturing the GT as a freezing in an 
aperiodic crystal structure. In this approach system relaxation 
towards equilibrium can be thought as a series of activated 
jumps over energy barriers in its potential energy landscape^. 
Following RefJ 2 -^ we now give a derivation of an equation re- 
lating DW factor and structural relaxation, which is useful in 
the context of this paper. 

For the sake of simplicity we restrict to the one-dimensional 
case where two minima are separated by a distance 2a (see 
Fig. [TJ. Referring to Fig. [T| we expand the potential U (x) 
around the minimum on the left, whose position is labeled by 



U(x) = U + -(x-x 1 y 



(1) 
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Fig. 3 Self part of intermediate scattering function for different q 
from q m j„ to values around q max - Solid and dashed black lines are 
fits to a stretched exponential with p = 0.58 (dashed) and (3 = 1.4 
(solid) for minimum and maximum wave vectors considered here 
respectively. 

Since system relaxation requires getting over the energy bar- 
rier AE, if T a is the system relaxation time and To the micro- 



scopic time: 



AE 

T« = T exp — 
k R T 



From Eq. ([TJ we can express AE as: 

AE = -a 2 
2 

and from equipartition theorem: 

k B T = A(u 2 ) 



(2) 



(3) 



(4) 



where (u 2 ) = (x 2 ). Inserting Eqs. Q and ^ into Eq. |2ji one 
obtains finally: 

T« = T exp ( — ) (5) 



It is important to note that Eq. <j5j is expected to fail if the 
amplitude of rattling motion (u 2 ) becomes comparable to a 2 . 

A natural generalization of Eq.([5]) can be achieved adopt- 
ing a proper distribution p(a 2 ) of the squared displacement a 2 
needed to overcome energy barriers, i.e. in our present study 
to break bonds. We note that the squared displacement a 2 
is the cumulative displacement of the particles that move^H 
hence according to Central Limit Theorem a suitable choice 
for p{a 2 ) is a truncated gaussian form, i.e. 



p(a 2 ) = 



A exp 



(a 2 -a 2 ) 2 
2a 2 



if a > a,„i„ 
otherwise 



(6) 



where A is a normalization factor and a 2 nin is the minimum 
displacement to reach the transition state. Averaging the Eq. 
|5]) over the distribution given by Eq. (|6j, the following gen- 
eralized HW equation is obtained: 
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(7) 



The gaussian form for p(a 2 ) is supported also by other con- 
siderations. For example if we substitute back kgT into Eq. 
(j7]) using Eq. Q we end up with the following equation: 



T a = To exp 



A 2 a 2 



Aa 2 



2k B T %{k B T) 2 



(8) 



Experimental data for both supercooled liquids^] and poly- 
mers iazj together with theoretical approaches support gaus- 
sian form for p(a 2 ) in Eq. <JsJ . Furthermore using Eq. |3j 
into Eq. ^ to eliminate a 2 a gaussian distribution for energy 
barriers is attained in accordance with other studies^!. 

These ideas have been originally developed for the glassy 
dynamics of dense systems, where the caging occurs due to 
the high density. In the following we would like to try and 
apply them to the glassy dynamics of colloidal gels, where 
some hints of caging phenomena appear, although densities 
can be very low. The caging in these systems rather originates 
from the formation of persistent bonds 33 , which eventually 
lead to an interconnected network structure. 

3 Methods 
3.1 Model 

We refer to the studies carried on in RefsJ^SHSl the colloidal 
gel model considers identical particles of diameter o interact- 
ing via a phenomenological potential V e ff, designed to ac- 
count for the presence of directional interactions. In gelling 
colloidal suspensions there are in fact several possible sources 
of anisotropic effective interactions, since the particle sur- 
face may not be smooth or the building blocks of the gel 
are not the primary particles but larger aggregates of irreg- 
ular shaped Confocal microscopy images obtained in re- 
cent experiments^^! confirm this scenario: the distribution 
of the particle coordination number n in very diluted gel net- 
works is strongly peaked around n ~ 2, 3. In the chosen model, 
therefore, the interaction potential is given as the sum of three 
different contributions, V e ff = Vzj +V4 + V3, where Vlj is a 
Lennard- Jones type of potential producing a narrow attrac- 
tive well, and Vd + V3 contains directionality and rigidity of 
inter-particle bonding 3 -^. Here we consider the same choice 



of parameters as in L ' 
investigated in^H 



and the range of volume fractions as 



3 



As reported in the previous studies, in this model at low 
temperatures the system aggregates into an open persistent 
network of chains connected by a few bridging points (nodes). 
This takes place via a random percolation mechanism, but 
once a percolating structure is formed, it rapidly evolves to- 
wards a persistent, fully connected open network. The forma- 
tion of the persistent network produces the coexistence, in the 
gel, of very different relaxation processes at different length 
scales: the relaxation at high wave vectors is due to the fast co- 
operative motion of pieces of the gel structure (e.g. the chains 
connecting two nodes), whereas at low wave vectors the over- 
all rearrangements of the heterogeneous gel make the system 
relax via a stretched exponential decay of the time correla- 
tors. The coexistence of such diverse relaxation mechanisms 
is characterized by a typical crossover length which is of the 
order of the network mesh size. The slow glassy dynamics at 
low wave vectors results to be directly connected to the pres- 
ence of cooperative processes which can be recognized, for 
example, in the rearrangements of the network nodes along 
the complex structure of the network itself 38 . This scenario is 
in agreement with the results of other recent studies on model 
colloidal gels^El. 



3.2 Simulations details 



We have used a MD code where the potential V e ff has been 
implemented via a suitable combination of the algorithms 
RATTLE and SHAKFM The unit of time is ^mo 2 /e, with 
m the mass of a particle and the data reported here refer to a 
time step of 0.002. The data refers to micro-canonical simu- 
lations performed with 8000 particles in cubic boxes of size 
L = 37.64,43.09,55.10 in unit of o, corresponding respec- 
tively to particle densities of p = 0.15, 0.1, and 0.05, i.e. to 
approximately volume fractions (]) ~ 0.075, 0.05, and 0.025. 
In the simulations 5 to 8 independent samples have been equi- 
librated starting from initial high temperature random config- 
urations by replacing particle velocities with values extracted 
from a Maxwell-Boltzmann distribution every A time steps 
(where A varied with temperature from 10 to 10 3 MD steps). 
After equilibration the energy is constant, showing no signifi- 
cant drift over the simulation time window, and different one- 
and two- time autocorrelation functions display the equilib- 
rium behavior, i.e. do not show any sign of aging. The data 
production starts from these equilibrated samples: the equili- 
bration time grows with the relaxation time in the system and 
at the lowest temperatures equilibration required up to 2 • 10 7 
MD steps. 



4 Results and Discussion 

4.1 Relaxation and Transport Properties 

We use the static structure factor S(q), defined as follows: 

^)^L< e ' q(r '~ r ')> (9) 

U 

to quantify the extent of spatial correlation in the system and 
obtain informations on the gel structure. In Fig. [2] (from the 
data of Ref. 38 ) S(q) of the gel network (i.e. at the lowest tem- 
perature considered) displays a peak around q max ss 8 corre- 
sponding roughly to the particle diameter. This peak basically 
arises from excluded volume interactions between particles, 
i.e. it approximately corresponds to the first peak of radial 
distribution function. In glassy systems the slow relaxation 
arises first, and has its strongest signature, at these wave vec- 
tors. It is clear from the figure that in the gel significant spatial 
correlations are present also at smaller wave vectors. In Refill 
S(q) has been compared to the static structure factor of a poly- 
mer chain solutionis! length scales matching smaller wave 
vectors 2.0 < q < 7.0 can interpreted as an intra-molecular 
regime for spatial correlations of the aggregates (i.e. chains). 
Mesoscopic and macroscopic length scales q < 2.0 can be in- 
stead thought of as corresponding to inter-molecular regime, 
due to the long-range interactions induced by the formation 
of the persistent gel network. Correlations in the particle mo- 
tion over different length scales can be effectively quantified 
in terms of ISF: 

F MS) = ^t(e'^-^) (10) 
The analysis of its 

behaviorism 

indicates that in the gel the 
slowest modes correspond to the inter-molecular regime of 
wave vectors. In Fig. ^/^(g,?) is plotted as a function of the 
time, rescaled by the relaxation time X a (q), at the lowest tem- 
perature T — 0.05 for different wave vectors. X a (q) has been 
calculated from F s (q,x a (q)) = 1 /<?. The figure well shows that 
the stretched exponential decay exp{— [f/T a (g)]P} with p < 1. 
typical of glassy dynamics, arises only at low wave vectors 
(q<\ .0). In other words, the slow structural relaxation modes 
of the gel structure can be detected only at low q. In contrast, 
at higher wave vectors, the time decay of correlations is faster 
than exponential (P ~ 1.4): in Refs J 37 ! 63 ! these processes have 
been carefully analyzed and it has been shown that they are 
due to fast coherent motion of pieces of the gel network (i.e. 
the chains between two bridging point or nodes). 

The overall scenario of relaxation modes in the ge l is there- 
fore rather different from the one discussed irJ 25 l 43 l for glassy 
systems at high densities. It is interesting to notice that here 
the particle bonding is the basic ingredient of the structural ar- 
rest, since the persistence of the gel network certainly relies 
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upon the long living bonds. On the other hand, the analysis 
of the relaxation modes well elucidate that formation of sin- 
gle bonds cannot be responsible, on its own, for the coopera- 
tive glassy dynamics which instead arises from the long-range 
correlations between them induced by the network^. 

Within this picture, the question of the existence and na- 
ture of a direct correlation between the single particle average 
localization and the structural arrest, as discussed in Sec |4.2| 
becomes particularly intriguing and is the main focus of this 
work. 

To this aim, we have calculated the relaxation time T a asso- 
ciated to the structural relaxation of the gel as x a = T a (q m i„), 
where q m j n is the smallest wave vector compatible with pe- 
riodic boundary conditions in our simulations boxes, i.e. it 
corresponds to length scales of the order of the box size. 

For evaluating the rattling amplitude of the particle motion 
corresponding to caging, we consider the MSD: 



(^(0> = ~£«)-ri(o)||; 



(11) 



Figure [4] shows MSD divided by time for all temperatures at 
volume fraction = 0.05 (from Ref.^j). Since the system is 
very diluted, localization phenomena can be very weak and 
this type of plot helps to better recognize different regimes. 
At very short times a ballistic regime is found where MSD in- 
creases according to (r 2 (t)) ~ (3k[,T /m)t 2 (i.e.(r 2 (f))/f °< t). 
Formation of bonds with other particles slows down the dis- 
placement and (r 2 (t))/t shows an inflection point. At high 
temperatures, bonds break within a time interval much smaller 
than T a and particle starts diffusing, i.e. (r 2 (t))/t eventually 
reaches a plateau. At T < 0. 1 bond lifetime sets instead the 
longest relaxation time-scale in the system and in this regime 
the formation of the gel network starts, with the MSD becom- 
ing increasingly sub-diffusive over times much longer than the 
localization process related to the rattling of the particle within 
the bonding length scale. Therefore we evaluate the caging 
from this first localization process as explained in the follow- 
ing. 

4.2 Scaling between relaxation and caged dynamics 

Following the discussion in RefJ^S we evaluate the DW factor 
in our gel system in order to provide a suitable characteristic 
length scale for the particle temporarily trapped into the cage 
due to bonds formation. DW factor can be defined picking 
a suitable value of MSD within a time window that begins 
just after ballistic regime and that ends just before structural 
relaxation sets in. First we have to identify such time window 
and to do that we consider the slope A(f ) of MSD in a log-log 
plot, i.e.: 

3log(r 2 ) 



Af 



9 log* 



(12) 



Representative plots of A(f ) for our gel system can be found in 
the inset of Fig. [4] The short-time ballistic regime corresponds 
to A(f ) « 2 while the long time diffusive regime corresponds to 
A(t) s» 1. Between these two regimes an intermediate regime 
is present where caging of particles due to bonds gives rise to 
a clear minimum of A(f ) (see Fig. |4}. We thus define the DW 
factor (u 2 ) as follows: 



(u 2 ) = (r 2 (t = n) 



(13) 



where t* is the time corresponding to the minimum of A(f) 
within this intermediate regime. 
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Fig. 4 (r 2 (t)} ft as a function of time at (j) = 0.05. Inset: Logarithmic 
derivative showing the minimum corresponding to caging. The 
arrows indicate the regime chosen for the evaluation of DW. 

We are now in a position to establish a correlation between 
structural relaxation and caged dynamics. In Fig. |5jlog(x a ) is 
plotted against the inverse of DW factor 1 / (u 2 ) for the three 
volume fractions investigated. The figure clearly shows that 
X a is strongly correlated to the DW factor. The correlation 
has also a specific form (see the fitting curve in the figure) 
well agreeing with the prediction of Eq. This is far 

from being obvious because here the structural relaxation is 
related to length scales of the order of simulation box (i.e. 
2n/q m i n ) whereas the DW factor corresponds to caging phe- 
nomena occurring on much smaller length scales of the order 
of first neighbor distance, i.e. q ss 2%/a « q max (see Fig. |2j), 
where the relaxation is rather dominated by other mechanisms 
(see Fig(3]). It is also remarkable that the all data for differ- 
ent volume fractions collapse onto the same scaling curve, 
because the structure, as described by spatial correlations in 
Fig|2j changes significantly with at small q (S(q m i„) in- 
creases in fact by almost an order of magnitude from = 0.075 
to = 0.025), whereas it is not affected by changes in (j) at 
large wave vectors q. 
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Fig. 5 Scaling of gel data for all densities investigated, i.e. 

= 0.025, 0.05, 0.075. Dashed line is a fit of all data to the function 

logT a = a+P(« 2 )- 1 +j(u 2 )- 2 with a = 1.893, P = 0.0177 and 

7=0.00144. 



Fig. 6 Plot of log(l a {q) ■ q lA ) versus 1 /u 2 for different q at 
(]) = 0.075. Upon increasing the wave vector from q m i„, the data 
strongly depart from the scaling form of Eq[7] 



where: 



The correlation between DW factor and structural relax- 
ation time T a points to a strong correlation between the long- 
range correlations established by the network and the localiza- 
tion within the structure at the level of the single particle. 

To better elucidate the nature of the scaling found, we 
have also investigated the existence of correlations at differ- 
ent wave vectors. In the gel network, T a ( q) dis plays a com- 
plex dependence on q as discussed in Refs J 37 l 38 l In particular, 
T a (q) °c q~ 1A at large q, corresponding to the regime where 
relaxation is dominated by fast collective motion of pieces of 
the structures (i.e. chains between two nodes). Therefore in 
Fig|6]we have used the same type of plot as in Fig |5] where 
T a (<?) has been rescaled with q 1A . The data refer to different q 
at (j) = 0.075. At the largest q the data well collapse on top of 
each other and display a very different dependence on 1/ (u 2 ). 
Upon decreasing q the data depart from the q 1 scaling in q 
and also approach the scaling form of Fig|5] This analysis fur- 
ther demonstrate that the scaling found specifically relates the 
average localization at the level of the single particle (as quan- 
tified by (u 2 )) to the relaxation modes associated to the arising 
of glassy, cooperative dynamics. 

4.2.1 Comparison with other MD Studiesln Ref.El it 
has been shown that for several model glassformers and ex- 
perimental systems plotting log(x a ) versus (u l)/{u 2 ) (where 
(u 2 ) is the DW factor at the GT) all data scale onto the same 
master curve, i.e.: 



■ (u 
logT„ = a + P 




a 
P 



-0.424(1) 
? 

21nl0(«2) " 

a 2 , 
81nl0( M 2)2 



1.62(6) 



12.3(1) 



(15) 
(16) 

(17) 



(14) 



The scaling form obtained in Fig. [5] for the colloidal gel 
can in fact be superimposed on the universal curve of Eq. ( [Pi} 
with a suitable vertical shift (a' = a — 2.33) and upon using on 
jc-axis the scaled variable (u 2 }/(u 2 ) with (u 2 ) = 0.104 in order 
to rule out any trivial dependence on time and length scales. 
Using such vertical shift and such value for (u 2 ) to adjust gel 
data we can compare them to the results obtained from soft 
binary mixture sE3 and polymer systems 25 as shown in Fig j7] 
It is clear from this figure that within the accuracy (marked 
by solid lines) the scaling procedure works well also for the 
colloidal gel model considered in the present paper. 



5 Conclusions 

We have investigated possible correlations between the local- 
ization at the level of the single particle and structural relax- 
ation in a model colloidal gel at very low volume fractions, 
with directional effective interactions and local rigidity. We 
have found that strong correlations are present over different 
length scales. In this type of systems, the localization of par- 
ticles due to persistent bonding is of course the initiator of the 
process that leads to gelation, but cannot be responsible, on its 
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Fig. 7 Scaling of the structural relaxation time vs the reduced DW 
factor of polymers 25 , soft binary mixtures (BM) 43 and colloidal gel 
(present work). For the colloidal gel {u 2 ) l l 2 = 0.103 and data have 
been shifted vertically by —2.33 in log scale. Dashed line is the 
universal curve defined in Eq. \14) . Solid lines bound the accuracy 
of Eq [7] obtained fitting data from polymer models studied inland 
correspond to the two definitions (u 2 ) = (r 2 (t = 0.6)) (magenta) 
and (u 2 ) = (r 2 {t = 1.4)) (red) (see 43 for more details). 



own, for structural arrest, which arises thanks to the formation, 
eventually, of an interconnected network structure. Remark- 
ably, we have found that relaxation modes at the lowest wave 
vectors, i.e. over length scale much larger than the bond local- 
ization length, strongly correlate to the localization of single 
particles within the typical bonding length. We think that these 
findings indicate a type of feedback mechanism between dy- 
namical processes at different length-scales: particle bonding 
leads to the network formation and long range correlations in- 
duced by the presence of the network actually transforms the 
particle bonding into a glassy caging, coupling eventually the 
particle localization to the glassy structural arrest. Moreover, 
we have shown that the gel data display the same scaling form 
found for a large class of dense glassy systems in experiments 
and simulations (molecular glasses, polymers, etc. ): this re- 
sult strongly suggests that the scaling captures the essential, 
basic ingredients in the physics of glassy structural arrest and 
it is an extremely powerful tool for devising its possible uni- 
versal features. 
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